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Abstract 

We present an algorithm for the numerical calculation of one-loop QCD amplitudes. The 
algorithm consists of subtraction terms, approximating the soft, collinear and ultraviolet 
divergences of one-loop amplitudes and a method to deform the integration contour for the 
loop integration into the complex space. The algorithm is formulated at the amplitude level 
and does not rely on Feynman graphs. Therefore all required ingredients can be calculated 
efficiently using recurrence relations. The algorithm applies to massless partons as well as 
to massive partons. 
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1 Introduction 



Multi-jet final states play an important role for the experiments at the LHC. An accurate descrip- 
tion of jet physics is therefore desirable. Although jet observables can rather easily be modelled 
at leading order (LO) in perturbation theory, this description suffers several drawbacks. A lead- 
ing order calculation depends strongly on the renormalisation scale and can therefore give only 
an order-of-magnitude estimate on absolute rates. Secondly, at leading order a jet is modelled 
by a single parton. This is a very crude approximation and oversimplifies inter- and intra-jet 
correlations. The situation is improved by including higher order corrections in perturbation 
theory. 

At present, there are many next-to-leading order (NLO) calculations for 2 — > 2 processes at 
hadron colliders, but only a few for 3 or more partons in the final state [1-49]. It is desirable 
to have NLO calculations for 2 — > n processes in hadron-hadron collisions with n in the range 
of n = 3, ...,7. However, the complexity of the calculation increases with the number of final 
state particles. For any NLO calculation there are two parts to be calculated: the real and the 
virtual corrections. Almost without exceptions all examples cited above use the dipole formalism 
[50-53] to subtract out the infrared divergences from the real corrections. The subtracted real 
correction term is then integrable in four dimensions and can be calculated numerically by Monte 
Carlo techniques. By now there are several implementations for the automated construction 
of subtraction terms [54-59]. The required Born amplitudes can be calculated efficiently with 
the help of recurrence relations [60-66]. The calculation of the virtual corrections for QCD 
processes with many external legs has been considered to be a bottle neck for a long time. The 
past years have witnessed significant progress in this direction. The main lines of investigation 
focus on a perfection of the traditional Feynman graph approach [67-77] or are based on unitary 
methods [78-92]. 

In this paper we would like to discuss a third and purely numerical approach. To this aim 
we extend the subtraction method to the loop integration of the virtual corrections and we eval- 
uate the subtracted virtual corrections numerically with the help of a suitable chosen contour 
deformation. The method is formulated in terms of amplitudes and does not rely on Feynman 
graphs. Therefore all ingredients can be calculated efficiently using recurrence relations. Purely 
numerical approaches have been discussed in the past [93-102]. The literature focuses either 
on individual Feynman graphs and subtraction terms for individual graphs or on a contour de- 
formation for infrared and ultraviolet finite amplitudes, where no subtraction terms are needed. 
Unfortunately the methods discussed in the literature for the subtraction terms on the one hand 
and for the contour deformation on the other hand are not compatible with each other and can- 
not be combined. Furthermore it is not clear if the methods discussed so far in the literature 
are sufficiently efficient to be applied to multi-parton processes. What is new in this paper is 
the development of compatible methods for the subtraction terms and the contour deformation 
and the combination of all relevant aspects into one formalism. Inspired by recent work on the 
structure of infrared singularities of multi-loop amplitudes [103, 104] we found that the soft and 
collinear subtraction terms can be formulated at the level of amplitudes, without referring to 
individual Feynman graphs [105]. This is a significant simplification and opens the door to an 
efficient implementation. In addition we need subtraction terms for the ultraviolet divergences. 
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In this paper we present a set of ultraviolet subtraction terms which have the expected form of 
local counterterms and which are particularly well suited for the numerical contour integration in 
the sense that this set does not introduce additional singularities along the contour. The second 
important ingredient of our method is an algorithm for the contour deformation. The subtraction 
terms eliminate only singularities, where the contour is pinched but leave singularities where 
a deformation into the complex plane is possible. The contour deformation takes care of these 
remaining singularities. It is highly non-trivial to find a general algorithm which avoids these sin- 
gularities and which leads to stable Monte Carlo results. We achieve this goal by first introducing 
Feynman parameters. After Feynman parametrisation the contour deformation of the loop mo- 
mentum is straightforward. In order to avoid all singularities we have to deform the integration 
over the loop momentum and the integration over the Feynman parameters. For an amplitude 
with a large number of external particles we take additional measures to improve the efficiency 
of the numerical Monte Carlo integration. Our method works for massless and massive particles. 

This paper is organised as follows: In section 2 we give an overview of the general ideas 
behind our approach. In a shortened form these ideas have been presented in [106]. In section 3 
we provide a complete list of all subtraction terms. The infrared subtraction terms have been 
given for the first time in [105], we list them here again with a minor modification. The minor 
modification in the collinear subtraction terms adapts the subtraction terms to the chosen method 
of contour deformation. The ultraviolet subtraction terms are new. In section 4 we discuss in 
detail the contour deformation. Together with the subtraction terms this part constitutes the core 
of our method. In section 5 we discuss a few points which might help to understand our method 
or which might help to avoid possible pitfalls. Examples of a possible pitfall are diagrams like 
massless tadpoles, which give zero in an analytical calculation. These diagrams have to be in- 
cluded in the numerical calculation in order not to spoil the local cancellation of singularities. 
Section 6 discusses checks and simple examples. NLO results on more complicated processes 
will be published in a separate publication. Finally, section 7 contains a summary and our con- 
clusions. In an appendix we have documented certain technical aspects of our method. 

2 General setup 

In this section we give an overview of our method and define our notation. In subsection 2.1 
we start from the subtraction method for real corrections and present the extension to the virtual 
corrections. The pole structure in the dimensional regularisation parameter of one-loop QCD 
amplitudes is very well understood and recalled in subsection 2.2. Throughout this paper we 
work with colour ordered amplitudes. These are defined in subsection 2.3. Subsection 2.4 intro- 
duces the notation which we use for the kinematics. Our method works at the level of amplitudes. 
These can be calculated efficiently with the help of recurrence relations without relying on Feyn- 
man graphs. Recurrence relations are discussed in subsection 2.5. For the construction of the 
subtraction terms we need to know from which integration region divergences arise. This is 
reviewed in subsection 2.6. 
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2.1 The subtraction method 



The starting point for the calculation of an infrared safe observable O in hadron-hadron collisions 
is the following formula: 

1 



(O) = £ / dx\f a {x\) / dx 1 f h {x 1 )—— — - 

^j,J J 2^(5)n S pi n (l)n S pi n (2)/7 C olour(l)n C olour(2) 

/ d§ n (pi,p2;p3,---,Pn+2)0(pi,...,p n+ 2)\A„+2\ 2 ■ (1) 

n 

In this equation we have written explicitly the sum over the flavours a and b of the two partons 
in the initial state. In addition there is a sum over the flavours of all final state particles, which is 
not shown explicitly. The momenta of the two incoming particles are labelled p\ and p2, while 
P3 to Pn+2 denote the momenta of the final state particles. f a (x) gives the probability of finding 
a parton a with momentum fraction x inside the parent hadron h. 2K(s) is the flux factor, for 
massless partons it is given by 2K(s) = 2s. The quantity w S pin(0 denotes the number of spin 
degrees of freedom of the parton i and equals two for quarks and gluons. Correspondingly, 
n co i our (i) denotes the number of colour degrees of freedom of the parton i. For quarks, this 
number equals three, while for gluons we have eight colour degrees of freedom. The matrix 
element \^ n +2^ is summed over all colours and spins. Dividing by the appropriate number 
of degrees of freedom in the initial state corresponds to an averaging. d§ n is the phase space 
measure for n final state particles, including (if appropriate) the identical particle factors. The 
matrix element \& n +2\ 2 is calculated perturbatively. 

The contributions at leading and next-to-leading order are written as 



(0} L0 = fo n do B , 



(0) NLU = J O n+i do K + J O n da v + J O n do L . (2) 

n+l n n 

Here a rather condensed notation is used. do B denotes the Born contribution, whose matrix 
elements are given by the square of the Born amplitudes with (n + 2) partons l-^^l 2 - Similar, 
do R denotes the real emission contribution, whose matrix elements are given by the square of the 
Born amplitudes with (n + 3) partons l-^^l 2 . d<5 V gives the virtual contribution, whose matrix 
elements are given by the interference term of the one-loop amplitude ^ + \, with (n + 2) partons, 

with the corresponding Born amplitude ^^ 2 - denotes a collinear subtraction term, which 
subtracts the initial state collinear singularities. Taken separately, the individual contributions at 
next-to-leading order are divergent and only their sum is finite. In order to render the individual 
contributions finite, such that the phase space integrations can be performed by Monte Carlo 
methods, one adds and subtracts a suitable chosen piece [50-53]: 



[0) NL0 = j (On+idaK-Ondo^j + J [ 0„d<5 V + O n do c + O n j do* J . 



(3) 



n+l 
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The term (O n +\do R — 0„do A ) in the first bracket is by construction integrable over the (n + 1)- 
particle phase space and can be evaluated numerically. The subtraction term can be integrated 
analytically over the unresolved one-particle phase space. Due to this integration all spin- 
correlations average out, but colour correlations still remain. In a compact notation the result 
of this integration is often written as 

do c + Jdo A = I®do B + K®d<5 B + P®d<5 B . (4) 
i 

The notation ® indicates that colour correlations due to the colour charge operators T ; still re- 
main. The action of a colour charge operator T ; for a quark, gluon and antiquark in the final state 
is given by 

quark: **{..&...) (Tfj) * (...qj...) , 
gluon: sl*(...g c ...)(if ab )si[...g b . 
antiquark : SX * (—qi—) (~Tf) Si (... qj..) . (5) 

The corresponding formulae for colour charge operators for a quark, gluon or antiquark in the 
initial state are 

quark : SI* (...<?/...) (-Tf) SI (...qj...) , 
gluon: R* (.../...) {if ab )*L (...«*...), 
antiquark: SI* (—qi—) (Tfj) SI (...qj..) . (6) 

In the amplitude an incoming quark is denoted as an outgoing antiquark and vice versa. The 
terms with the insertion operators K and P do not have any poles in the dimensional regularisation 
parameter and pose no problem for a numerical evaluation. The term I ® do B lives on the phase 
space of the n-parton configuration and has the appropriate singularity structure to cancel the 
infrared divergences coming from the one-loop amplitude. Therefore do v + 1 ® do 3 is infrared 
finite. We emphasise that this cancellation occurs after the loop integration has been performed 
analytically in D dimensions. do v is given by 

da v = 2Re (V 0) *.3 (1) ) o n dty n . (7) 
S\ W denotes the renormalised one-loop amplitude. It is related to the bare amplitude by 

Sl^j denotes the ultraviolet counterterm from renormalisation. The bare one-loop amplitude 
involves the loop integration 

= f^G W (9) 

bare J (2^)^ are ' 
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where denotes the integrand of the bare one-loop amplitude. In this paper we extend the 
subtraction method to the integration over the virtual particles circulating in the loop. To this aim 
we rewrite eq. (8) as 



(1) 
bare 



+ A 



1) 
CT 



A 



(1) 
bare 



A 



(1) 
soft 



A 



(!)_-(!) 



coll 



A 



uv 



(1)^(1)^(1)^(1) 



soft 



l coll 



L UV 



(10) 



The subtraction terms A^l, ^^011 an< ^ -^uv are cnosen sucn that they match locally the singular 
behaviour of the integrand of #viL m ^ dimensions. The first bracket in eq. (10) can therefore 
be integrated numerically in four dimensions. The term A ^2 approximates the soft singulari- 
ties, A^l a PP rox i ma tes the collinear singularities and the term A^y approximates the ultraviolet 
singularities. These subtraction terms have a local form similar to eq. (9): 



A 



(1) 
soft 



d D k (i) 
(2ji) ZJ ^ soft ' 



j3 



(i) 

coll 



<nt (i) 

(27Tp^ co11 ' 



UV 



(my 



The contribution from the terms in the first bracket of eq. (10) can be written as 



/ 



2 Re 



A 



(or 



A 



bare 



A 



soft 



A 



coll 



A 



(1) 

uv 



dty n 



r 2Re 



A 



r (l) 

"bare 



1) 
soft 



1) 
coll 



^UV 



o„ + (e) . 



(12) 



The integral on the right-hand side is finite. It is one of the key ingredients of the method 
proposed here, that this rather complicated and process-dependent integral can be performed nu- 
merically with Monte Carlo techniques. We recall that the error of a Monte Carlo integration 
depends on the variance of the integrand and scales with the number of integrand evaluations 
N like 1 / v^V- It is important to note that the error does not depend on the dimension of the 
integration region. Eq. (12) gives for an observable a contribution to the next-to-leading order 
prediction. The right-hand side corresponds to a (3m) -dimensional integral. (The phase-space in- 
tegral is (3m — 4) -dimensional, the loop integral 4-dimensional.) In practise this (3m) -dimensional 
integral is done with a single Monte Carlo integration. There is no need to evaluate for a given 
phase-space point the inner four-dimensional loop integral by a separate Monte Carlo integration. 
This is essential for the efficiency of the method. 

The building blocks of the subtraction terms are process-independent. When adding them 
back, we integrate analytically over the loop momentum k. The result can be written as 



2 Re 



A 



(0) : 



^CT ^ -Soft ^ -*coll ^ ^UV 



(13) 



The insertion operator L contains the explicit poles in the dimensional regularisation parameter 
related to the infrared singularities of the one-loop amplitude. These poles cancel when combined 
with the insertion operator I: 



(I-f-L) ®do L 



finite. 



(14) 



The operator L contains, as does the operator I, colour correlations due to soft gluons. 
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2.2 Explicit poles in the dimensional regularisation parameter 

The explicit poles in the dimensional regularisation parameter £ of the individual pieces are well 
known. These poles are either of ultraviolet or infrared origin. Let us consider an amplitude 
with n q external quarks, nq external anti-quarks and n g external gluons in massless QCD. We set 
n = riq + riq +n g . Obviously we have n q = nq. After an analytical integration the poles in the 
dimensional regularisation parameter £ of a bare massless one-loop QCD amplitude are given by 



(i) 

bare 



47tr(l-e) 



2) Po 



I 2 T 2 £ 



-ZPiPj 

V 2 



^ 0) + O(£°). 



(15) 



Po is the first coefficient of the QCD |3-function and given by 

Po = —C A --T R N f . 



(16) 



The constants Ji are given by 



3^ L 



(17) 



The colour factors are as usual 

C A =N C1 C f 



2N r 



Tr 



1 

2' 



(18) 



The poles of eq. (15) are cancelled by and the infrared poles obtained from integrating the 
real emission contribution over the unresolved phase space. The leading order QCD amplitude 

/ 2") /2 

with n partons is proportional to <xy . For a massless QCD amplitude with n partons the 
ultraviolet counterterm is given by 



(l) 

CT 



a s (n - 2) (3q (q) 



4% 



(19) 



The insertion operator I contains the infrared poles obtained from integrating the real emission 
contribution over the unresolved phase space. The insertion operator I is given in massless QCD 
by 



a* 



2jtr(l-e) 



1 n 




+ o(e). (20) 



with 



K n — K„ 



18 6 ) A 9 f 



(21) 
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2.3 Colour decomposition 

Amplitudes in QCD may be decomposed into group-theoretical factors (carrying the colour 
structures) multiplied by kinematic functions called partial amplitudes [107-111]. These partial 
amplitudes do not contain any colour information and are gauge invariant objects. The colour 
decomposition is obtained by replacing the structure constants f abc by 



if abc = 2 



Tr ^f ar p^ r p c ^j Tr ^jbjajt 



(22) 



which follows from [T c \ T b ] = if abc T c . In this paper we use the normalisation 

TrT a T b = }_ d ab (23) 

2 

for the colour matrices. The resulting traces and strings of colour matrices can be further simpli- 
fied with the help of the Fierz identity : 

TfjTfi = Ifajk-^iA^j. (24) 

There are several possible choices for a basis in colour space. A convenient choice is the colour- 
flow basis [54, 1 12, 1 13]. This choice is obtained by attaching a factor 

V2T? (25) 

to each external gluon. As an example we consider the colour decomposition of the pure gluon 
amplitude with n external gluons. The colour decomposition of the tree level amplitude may be 
written in the form 

(\ n—2 
-jz) E 8 iai^a Sidles •■■ 8 <oi,M A » °(S°i'-> (26) 
v2/ ceS „/z„ 

where the sum is over all non-cyclic permutations of the external gluon legs. The quantities 
A n (g ai , ...,g 0n ), called the partial amplitudes, contain the kinematical information. They are 
colour ordered, e.g. only diagrams with a particular cyclic ordering of the gluons contribute. For 
the convenience of the reader we have listed the colour ordered Feynman rules in appendix A. 
Similar decompositions exist for all other Born QCD amplitudes. As a further example we give 
the colour decomposition for a tree amplitude with a pair of quarks: 

£ \h l Kjo 2 ---K n _ 2 j^n\q,go 1 ,---,g<5 n ^,q), (27) 



where the sum is over all permutations of the external gluon legs. In squaring these amplitudes a 
colour projector 



has to be applied to each gluon. At the one-loop level the colour decomposition is slightly 
more involved. Let us start with an example: The colour decomposition of the one-loop gluon 
amplitude into partial amplitudes reads [111]: 

An\gl,82,-,gn) = ( A= ) { £ N d^ j 8; _/ ...S^y A$ (# 0l , ...,g a „) 

VV2/ [aeS n /z n 

+ L L [^j^-KmjaJ {\ n+ ijc m+2 ---^c n ja m+l 

m oeS„/(Z m xZ„- m ) 

An)m {g(3\ i •■•i8a m >8c m+ i j ■■■i8a n ) | • (29) 

The one-loop partial amplitudes A„)n are again gauge invariant. The subleading partial ampli- 
tudes are related to the leading ones by 

A nM8h-,8m,8m+h-,8n) = (-1)'" £ A l|o > ->£a«) ' (30) 

C€{m,...,l}UI{m+l,...,n} 

where the sum is over all shuffles of the set {m, ...,1} with the set {m+ 1, It is therefore 

sufficient to focus on the leading partial amplitudes A^q. The leading partial amplitudes can be 
decomposed further into smaller objects called primitive amplitudes. The primitive amplitudes 
are separately gauge invariant. For the example of the one-loop gluon amplitude we have the 
decomposition of the leading partial amplitude into two primitive amplitudes, which are charac- 
terised by the particle content circulating in the loop: 

For the primitive amplitude A^q 1c there is either a gluon or a ghost circulating in the loop, while 

for the primitive amplitude A^l nf there is a quark circulating in the loop. 

Similar decompositions exist for all other one-loop QCD amplitudes. If the external legs of a 
one-loop QCD amplitude involve a quark-antiquark pair we can distinguish the two cases where 
the loop lies to right or to the left of the fermion line if we follow the fermion line in the direction 
of the flow of the fermion number. We call a fermion line "left-moving" if, following the arrow 
of the fermion line, the loop is to the right. Analogously, we call a fermion line "right-moving" 
if, following the arrow of the fermion line, the loop is to the left. Examples of diagrams regarding 
these types of primitive amplitudes are shown in fig. (1). It turns out that in the decomposition 
into primitive amplitudes a specific quark line is, in all diagrams which contribute to a specific 
primitive amplitude, either always left-moving or always right-moving [114]. Therefore in the 
presence of external fermions primitive amplitudes are in addition characterised by the routing of 
the fermion lines through the amplitude. In summary we can always write a full one-loop QCD 
amplitude as a linear combination of primitive amplitudes: 

*£> = IC,A« t . (32) 

j,k 
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Figure 1: Examples of diagrams regarding left-moving and right-moving primitive amplitudes: 
Diagrams (a) contributes to the left-moving primitive amplitude, while diagram (b) contributes 
to the right-moving amplitude. 



The colour structures are denoted by Cj, while the primitive amplitudes are denoted by A n j k . In 
the colour-flow basis the colour structures are linear combinations of monomials in Kronecker 
8 ; /s. In order to construct a primitive one-loop amplitude one starts to draw all possible planar 
one-loop diagrams with a fixed cyclic ordering of the external legs, subject to the constraint that 
each fermion line is either only left-moving or only right-moving. This set of diagrams is then 
further divided into the subset of diagrams with a closed fermion line and the ones without. The 
set of diagrams with a closed fermion line form a separate primitive amplitude. 

In the following we will work exclusively with primitive amplitudes. In order to simplify the 
notation we will drop the subscripts and simply write 

A (1) (33) 

for a primitive one-loop amplitude. The full one-loop amplitude is just the sum of several primi- 
tive amplitudes multiplied by colour structures. We stress a few important properties of primitive 
one-loop amplitudes: 

1 . Primitive amplitudes are gauge invariant. 

2. Primitive amplitudes have a fixed cyclic ordering of the external legs and a definite routing 
of the external fermion lines. 

3. The flavour of each propagator in the loop is unique: Either it is a quark propagator or a 
gluon/ghost propagator. 

Our method exploits these facts. The first property of gauge invariance is crucial for the proof of 
the method. The second point ensures that there are at maximum n different loop propagators in 
the problem, where n is the number of external legs. 

2.4 Kinematics 

We introduce some notation which will be used throughout this paper. Let be a primitive 
amplitude with n external legs. Since the cyclic ordering of the external partons is fixed, there 
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k n —i 



Pi Pn 

Figure 2: The labelling of the momenta for a primitive one-loop amplitude. The arrows denote 
the momentum flow. 

are only n different propagators occurring in the loop integral. With the notation as in fig. (2) we 
define 

j 

kj = k-qj, qj = Y,Pl- ( 34 ) 
1=1 

We can write the bare primitive one-loop amplitude as 

A ^ e = !iw G ^ G ^ =m %q-l? j +ib' (35) 

P(k) is a polynomial in the loop momentum k. The +z'8-prescription in the propagators indicates 
into which direction the poles of the propagators should be avoided. We denote by I g the set of 
indices j, for which the propagator j in the loop corresponds to a gluon. If we take the subset of 
diagrams which have the gluon loop propagator j and if we remove from each diagram of this 
subset the loop propagator j we obtain a set of tree diagrams. After removing multiple copies of 

identical diagrams this set forms a Born partial amplitude which we denote by 

We introduce two matrices, which depend on the external momenta (and the internal masses), 
but not on the loop momentum k. The kinematical matrix S is a n x n-matrix and is defined by 

Sij = (qi-qj) 2 -m- -mj. (36) 
The Gram matrix is a (n — 1) x (n — 1) -matrix and is defined by 

G\f = 2 (?i - 9a) ■ (qj ~ q a ) , (37) 
where the indices i and j take the values i,j^a. 
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2.5 Recurrence relations 



Although we use in the proof of our method the fact that a primitive amplitude can be written 
as a sum of Feynman diagrams, it is important to note that in the final formulae, which enter the 
numerical Monte Carlo program, only amplitudes occur - and not individual Feynman diagrams. 
There are several possibilities how these amplitudes can be calculated. A particular efficient 
method is based on recurrence relations. We first review the recursive method for tree-level 
partial amplitudes and discuss afterwards the necessary modifications for the computation of the 
integrand of a one-loop primitive amplitude. 

We start with the computation of tree-level partial amplitudes. Berends-Giele type recurrence 
relations [60] build tree-level partial amplitudes from smaller building blocks, usually called 
colour ordered off-shell currents. Off-shell currents are objects with n on-shell legs and one 
additional leg off-shell. Momentum conservation is satisfied. It should be noted that off-shell 
currents are not gauge invariant objects. Recurrence relations relate off-shell currents with n legs 
to off-shell currents with fewer legs. As an example we discuss the pure gluon current, which 
can be used to calculate the pure gluon amplitude. The recursion starts with n = 1 : 

/V) = ej(p). (38) 

8j* is the polarisation vector of the gluon corresponding to the polarisation X. The recursive 
relation states that in the pure gluon off- shell current a gluon couples to other gluons only via the 
three- or four-gluon vertices : 



P- 



+t E vr%(p] i ,..,p]%(P% + ;,..,P^Va(pt + ^..,pi") 

;'=i/=y+i 



,(39) 



where 



Pi J = Pi + Pi+i + .-.+Pj (40) 
and V3 and V4 are the colour ordered three-gluon and four-gluon vertices 



Vf pa = i(2^r-rg po -A vp ). (41) 
The recurrence relation is shown pictorially in fig. (3). The gluon current J p is conserved: 

ErfV* = °- < 42) 

From an off- shell current one easily recovers the on-shell amplitude by removing the extra prop- 
agator, taking the leg (n+ 1) on-shell and contracting with the appropriate polarisation vector. 
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Figure 3: The recurrence relation for the gluon current. In an off-shell current particle n + 1 is 
kept off-shell. This allows to express an off-shell current with n on-shell legs in terms of currents 
with fewer legs. 

Similar recurrence relations can be written down for the quark and antiquark currents, as well 
as the gluon currents in full QCD. The guiding principle is to follow the off-shell leg into the 
"blob", representing the sum of all diagrams, and to sum on the r.h.s of the recurrence relation 
over all vertices involving this off-shell leg and off-shell currents with less external legs. 

There are only a few modifications needed to compute the integrand of a one-loop primitive 
amplitude. We can divide all vertices into two classes: either a vertex is directly connected to 
the loop or it is not. Again we can follow the off-shell leg into the "blob". If the first vertex 
which we encounter belongs to the second class, then there is attached to this vertex a one- 
loop off-shell current with fewer legs. The other objects attached to this vertex are one or more 
tree-level off-shell currents. If on the other hand the first vertex belongs to the first class, then 
two edges of the vertex are connected to loop propagators. The object connected to these two 
edges is a tree-level off-shell current with two legs off-shell. It can be computed with methods 
similar to the computation for tree-level off-shell currents with one leg off-shell. We illustrate 
this principle in fig. (4) for the case of a toy theory with a single field and a single three-valent 
vertex. Within our method we also need an ultraviolet subtraction term. The complete ultraviolet 
subtraction term has a structure similar to ordinary ultraviolet counterterms: We can represent the 
ultraviolet subtraction term as a sum over diagrams, where each diagram has a tree- structure with 
exactly one propagator or vertex replaced by a basic ultraviolet subtraction term. The complete 
ultraviolet subtraction term can again be calculated recursively. As example we consider again 
the toy model with a single field and a single three-valent vertex. The recursion relation for the 
ultraviolet subtraction term is shown in fig. (5). 

2.6 Singular regions of the integrand of one-loop amplitudes 

In order to construct the subtraction terms we have to know in which regions of the integra- 
tion domain for the loop momentum divergences arise. There are three types of singularities in 
one-loop amplitudes, which have to be approximated by suitable local subtraction terms. The 
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Figure 4: The recurrence relation for the one-loop current in a toy model with a single field 
and a single three-valent vertex. The one loop currents are represented by an oval with a hole, 
tree-level currents are represented by an oval without a hole. 
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Figure 5: The recurrence relation for the ultraviolet subtraction terms in a toy model with a single 
field and a single three-valent vertex. Objects with a cross represent ultraviolet subtraction terms. 
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singularity types are soft, collinear and ultraviolet. Let us briefly review under which condi- 
tions these singularities occur [67,97, 115]. Soft singularities occur when a massless particle is 
exchanged between two on-shell particles. In the amplitude 

41 = 0*>n ( ^/ m , + , 

this corresponds to the case 

mj = 0, p] - m)_ x = 0, p 2 j+l - m 2 j+l = 0. (44) 

In that case the propagators (j — 1), j and (j' + l) are on-shell. The singularity comes from the 
integration region 

k ~ qj. (45) 

A collinear singularity occurs if a massless external on-shell particle is attached to two massless 
propagators. In the amplitude of eq. (43) this corresponds to 

pj = 0, mj-i = 0, mj = 0. (46) 

In that case the propagators (j — 1) and j are on-shell. The singularity comes from the integration 
region 

k ~ qj-xpj, (47) 

where x is a real variable. 

Finally ultraviolet singularities arise when components of the loop momentum tend to infin- 
ity. We will work in Feynman gauge throughout this paper. Therefore any loop integral with n 
propagators in the loop is maximally of rank n. Power counting arguments show immediately 
that all diagrams with five or more propagators in the loop are ultraviolet finite. Therefore all 
ultraviolet divergent diagrams have four or less propagators in the loop. It can be shown that 
the ultraviolet divergent diagrams are only those which are propagator or vertex corrections. Of 
course this has to be the case for a renormali sable theory. 



3 The subtraction terms 

In this section we present all subtraction terms for the numerical calculation of one-loop QCD 
amplitudes. In subsection 3.1 we give the infrared subtraction terms for massless QCD. The 
generalisation to massive QCD is presented in subsection 3.2. The ultraviolet subtraction terms 
for massless and massive particles can be found in subsection 3.3. All subtraction terms are 
added back in integrated form. The sum of all integrated subtraction terms (plus the ultraviolet 
counterterms) defines the insertion operator L. This is discussed in subsection 3.4. 
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Figure 6: The configuration for the soft limit. 

3.1 The infrared subtraction terms for massless QCD 

In this section we present the infrared subtraction terms for massless QCD. The extension of the 
infrared subtraction terms to massive particles is discussed in the next section. We start with the 
soft subtraction term, which is given by 

G soft - 1 L k 2 k 2 k 2 A j ■ < 48 ) 

jelg K j-i K j K j+i 

We recall that I g denotes the set of indices j, for which the propagator j in the loop corresponds 
to a gluon. 

The soft subtraction term is derived as follows: In the case where gluon i is soft, the corre- 
sponding propagator goes on-shell and we may replace in all Feynman diagrams which have the 
propagator i the metric tensor of this propagator by a polarisation sum and gauge terms: 

Here k\ denotes the on-shell limit of k; and denotes the sum over the physical polarisations: 

J •* 

d^(k,n) = ^(k^el^n) = -iT + 2 - . (50) 

rf is a light-like reference vector. We note that self-energy diagrams are not singular in the soft 
limit, therefore adding them to the loop diagrams will not change the soft limit. With the inclu- 
sion of the self-energy diagrams and a corresponding replacement as in eq. (49) the contribution 
from the polarisation sum in eq. (49) makes up a tree-level partial amplitude, where two gluons 
with momenta k^ and have been inserted between the external legs j and j +1. This is illus- 
trated in fig. (6). In the soft limit this tree-level partial amplitude is given by two eikonal factors 
times the tree-level partial amplitude without these two additional gluons: 
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In the soft limit we may replace 2pj ■ kj by kj_ { and 2pi+j ■ (—kj) by kj +l . Eq. (51) then leads 

to eq. (48). The terms with k^rC and rfk^ in eq. (49) vanish for the sum of all diagrams due to 
gauge invariance. Integrating the soft subtraction term we obtain 

r l2e f d°k (!) 1 y 2 ( -2p rPj+i y e (Q) 

We have multiplied the left-hand side by S^ l fi 2c , where S £ = (47l) e e _£Y£ is the typical volume 
factor of dimensional regularisation, Je is Euler's constant and /j is the renormalisation scale. 
Let us now consider collinear singularities. The collinear subtraction is given by 

= f i (~ 2 ) k y k 2 ' + k 4 ; Mf. <«) 

jei g \ K i i% K j K j+i 
where the symmetry factors are given by 

S q = S q =l, S 8 = ^. (54) 

The function guv ensures that the integration over the loop momentum of the subtraction terms 
in eq. (53) is ultraviolet finite. The function guv has the properties 

lim guv {kj- i,kj) = h lim guv {kj_ u kj) = O ( - ) . (55) 

There are many possible choices for the function guv- We use here a choice which is compatible 
with the contour deformation discussed in sect. (4). We take the function guv as 

k 2 _ k 2 

gw = 1-- \ l J - i2 . (56) 

Q is an arbitrary four- vector independent of the loop momentum k and is an arbitrary scale. 
Since these two quantities are arbitrary, there are no restrictions on them, they even may have 
complex values. We will later choose /j^y purely imaginary with Im /j uv < 0. This will en- 
sure that the denominator of eq. (56) does not introduce additional singularities for the contour 
integration. 

The collinear subtraction term is derived as follows: We have to consider configurations 
where two adjacent propagators go on-shell with a massless external leg in between. These 
configurations are shown in fig. (7). The diagrams (c) and (d) of fig. (7), where an external 
gluon splits into a ghost-antighost pair or into a quark-antiquark pair are in the collinear limit not 
singular enough to yield a divergence after integration. Therefore we are left with diagrams (a) 
and (b). Let us first consider the q qg splitting as in diagram (a). In Feynman gauge one can 
show that only the longitudinal polarisation of the gluon contributes to the collinear limit. The 
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Figure 7: Configurations for the collinear limit. Only diagrams (a) and (b) lead to a divergence 
after integration. Diagrams (c) and (d) are not singular enough to yield a divergence after inte- 
gration. 




Figure 8: In the collinear limit the sum of the singular diagrams in Feynman gauge of the ampu- 
tated one-loop amplitude is equal to the negative of the self-energy insertion on the external line 
where the gluon carries a longitudinal polarisation. 

same holds true for the g — > gg splitting of diagram (b). In this case the collinear limit receives 
contributions when one of the two gluons in the loop carries a longitudinal polarisation (but not 
both). The external gluon has of course physical transverse polarisation. It is well-known that 
the contraction of a longitudinal polarisation into a gauge invariant set of diagrams yields zero. 
The shaded "blobs" of picture (a) and (b) of fig. (7) consist almost of a gauge invariant set of 
diagrams. There is only one diagram missing, where the longitudinal polarised gluon couples 
directly to the other parton connected to the shaded blob. This is a self-energy insertion on an 
external line, which by definition is absent from the amputated one-loop amplitude. We can now 
turn the argument around and replace the sum of collinear singular diagrams by the negative of 
the self-energy insertion on the external line. This is indicated in fig. (8) for the q qg splitting 
and in fig (9) for the g — >■ gg splitting. The self-energy insertions on the external lines introduce a 
spurious 1 / /^-singularity. In order to calculate the singular part of the self-energies we regulate 
this spurious singularity by putting p 1 - slightly off-shell, but keeping kj_i and kj on-shell and 
imposing momentum conservation. We can use the same parametrisation as in the real emission 
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Figure 9: In the collinear limit the sum of the singular diagrams in Feynman gauge of the ampu- 
tated one-loop amplitude is equal to the negative of the self-energy insertion on the external line 
where one of the gluons carries a longitudinal polarisation. 



case: 



k | n 

kj_! = xp + k ± — -, 

x {2p-n) 

-kj = (l- x )p-k ± - k± H (57) 

(l-x){2p-n) 



with p L = n 2 = and 2p ■ k± = 2n ■ k± = 0. The singular part of the self-energies is proportional 
to 



« = -2k-^\-— X + 2 )' f - 

= _ * (A +2 ) (_f + 2 p^p1\ . (58) 



8 ^ 88 2kj- 1 -k j \ x 1-x J\ 2p-n 

The terms with 2/x and 2/(1 — x) correspond to soft singularities and have already been sub- 
tracted out with the soft subtraction term g[^ v In the collinear limit we therefore just have to 
subtract out the terms which are non-singular in the soft limit. These terms are independent of x 
and lead to eq. (53). If we compare the soft and collinear subtraction terms for the integrand of a 
one-loop amplitude with the subtraction terms for the real emission, we observe that there are no 
spin correlations in the subtraction terms for the integrand of the one-loop amplitude. In the real 
emission case spin correlations occur in the collinear limit. This can be understood as follows: 
From the proof of eq. (53) one can see that in the collinear limit always one of the collinear 
gluons carries an unphysical longitudinal polarisation. Hence, there are no correlations between 
two transverse polarisations of a gluon. 

Integrating the collinear subtraction term we obtain 



3.2 The infrared subtraction terms for massive QCD 

In this section we present the infrared subtraction terms for massive QCD. This is in particular 
relevant to top quark physics. There are only a few modifications necessary with respect to the 
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massless case, which will be discussed in the following. The modification for the unintegrated 
soft subtraction terms is straightforward: 



4 Pj-Pj+l 



,(0) 



m 



7+1 



(60) 



As before, the sum runs over all particles in the loop which are gluons. Integrating the soft 
subtraction term we have to distinguish whether the masses my_i and mj + i are zero or not. The 
result can be written as 



1 



* J (2n)D soft_ (47t)2r(l-e) y , 



Y,C({ P j + Pj+i) 2 ,m)_ x ,m) +l ,Li 2 )Af + o(e), (61) 
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(63) 



The modification for the collinear subtraction term is even simpler: There is no collinear singu- 
larity if an external quark or antiquark is massive. It suffices therefore to define Sq = Sq — for 
a massive quark or antiquark. 



3.3 The ultraviolet subtraction terms 

A primitive one-loop QCD amplitude contains, apart from infrared divergences, also ultraviolet 
divergences. In an analytical calculation regulated by dimensional regularisation these diver- 
gences manifest themselves in a single pole in the dimensional regularisation parameter £. If we 
would do the calculation within a cut-off regularisation we would find a logarithmic dependence 
on the cut-off. 

Within a numerical approach we need subtraction terms which approximate the ultraviolet 
behaviour of the integrand locally. We first note that the one-loop amplitude has, in a fixed 
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direction in loop momentum space, up to quadratic UV-divergences. These are reduced to a 
logarithmic UV-divergence only after angular integration. For a local subtraction term we have 
to match the quadratic, linear and logarithmic divergence. Since we work in Feynman gauge 
throughout this paper it follows that any loop integral with n propagators in the loop is maximally 
of rank n. Power counting arguments show immediately that all diagrams with fiver or more 
propagators in the loop are ultraviolet finite. Therefore all ultraviolet divergent diagrams have 
four or less propagators in the loop. It can be shown that the ultraviolet divergent diagrams are 
only those which are propagator or vertex corrections. The ultraviolet subtraction terms have to 
subtract out correctly the exact pointwise ultraviolet behaviour of the one-loop amplitude. There 
are several possibilities how the ultraviolet subtraction terms can be chosen. We present here a 
list of ultraviolet subtraction terms which have two additional properties: 

1 . The unintegrated ultraviolet subtraction terms have all singularities localised on a single 
surface 

(k-Qf-vlv = 0. (64) 

Q is an arbitrary four-vector independent of the loop momentum k. By choosing /jy V 
purely imaginary with Im < we can ensure that the contour for the integration over 
the loop momentum k never comes close to the singular surface defined by eq. (64). We 
remark that also the four- vector Q is allowed to have complex entries. We will choose Q 
as a function of the (real) external momenta and the (complex) Feynman parameters. 

2. We choose the ultraviolet subtraction terms such that the finite part of the integrated ul- 
traviolet subtraction terms is independent of Q and proportional to the pole part, with the 
same constant of proportionality for all ultraviolet subtraction terms. This ensures that the 
sum of all integrated UV subtraction terms is again proportional to a tree-level amplitude. 

To present the ultraviolet subtraction terms we follow the notation of fig. (2). This implies that we 
take all external momenta as outgoing and that the direction of the flow of the loop momentum is 
clock-wise. We label the external momenta clock-wise with pi, p2, etc.. With the conventions of 
fig. (2) the loop momentum k is the momentum of the loop propagator preceding the external leg 
with momentum p\. We set D = 4 — 2e for the number of space-time dimensions and we define 
khy 

k = k-Q. (65) 

It is sufficient to present the ultraviolet subtraction terms with these conventions. The general 
case can always be obtained by an appropriate substitution. Technically, the ultraviolet subtrac- 
tion terms are derived by expanding the loop propagators around the propagator corresponding 
to eq. (64). For a single propagator we have 

1 = 1 f | 2k-(p-Q) ( p -Q) 2 - m 2+^ y | [2k-(p-Q)} 2 \ 

+o (w) • (66) 
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(a) (b) (c) 

Figure 10: The master diagrams for the self-energies. Not shown are diagrams which are ob- 
tained from these by pinching some of the propagators. Also not shown are diagrams with a 
closed ghost loop. Diagram (a) gives the self-energy of the quark propagator, diagram (b) cor- 
responds to the leading colour self-energy of the gluon propagator while diagram (c) gives the 
^/■-contribution to the self-energy of the gluon propagator. 



We start with the ultraviolet subtraction term for the quark propagator. The relevant diagram 
whose large fc-behaviour has to be approximated is shown in fig. (10 a). For the quark propagator 
we take as UV- subtraction term 



-ffi, (1) 
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k-(p-iQ) y 



(67) 



The term proportional to /j^y in the numerator is not divergent, but ensures that the finite part of 
the integrated expression is proportional to the pole part. Integration yields 
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1 



uv 



(4jc) 
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+ 6 



(68) 



We see that the finite part is given by — ln^yy/ /j 2 ) times the pole part. With the choice /j^y = /j 2 
the finite part of the integrated ultraviolet subtraction terms will be zero. However this is not 
the way how we will use this formula. Our priority is to avoid additional singularities for the 
numerical loop integration. Therefore we will choose /j^y along the negative imaginary axis. On 
the other hand /u 2 is the usual renormalisation scale and conventionally taken as a positive real 
number. 

Next we consider the gluon propagator. We divide the ultraviolet subtraction terms into a 
leading colour part and a part proportional to Nf. An example of a diagram contributing to each 
part is shown in fig. (10 b) and fig. (10 c). It is convenient to define 



Rf 
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(69) 
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The leading colour UV- subtraction term is given by 
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Integration yields 
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(72) 



For the quark-gluon vertex there is a leading colour contribution and a subleading colour con- 
tribution. Representative diagrams are shown in fig. (11 a) and fig. (11 b). The subtraction term 
for the leading colour contribution is given by 
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.j 2e f d D k 
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A{\-t)kfW-2p 2 ]W f 
(k 2 -tiv) 3 



(73) 





Figure 11: The master diagrams for the vertex corrections. Not shown are diagrams which 
are obtained from these by pinching some of the propagators. Also not shown are diagrams 
with a closed ghost loop. Diagram (a) gives the leading colour vertex correction to the quark- 
gluon vertex, diagram (b) gives a subleading colour contribution. Diagram (c) corresponds to the 
leading colour vertex correction to the three-gluon vertex, diagram (d) gives the ^-contribution. 
Diagram (e) corresponds to the leading colour vertex correction to the four-gluon vertex, diagram 
(f) gives the A/y-contribution. The external momenta are always directed outwards, the loop 
momentum is always taken to flow clock-wise. 
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For the subleading colour contribution we have 
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Integration leads to 
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For the three-gluon vertex we have a leading colour contribution and a contribution proportional 
to Nf. Representative diagrams are shown in fig. (1 1 c) and fig. (lid). For the three-gluon vertex 
we define 
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The subtraction terms are given by 
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Integration leads to 

Cc = ^^fe P -^)+^ VP W-^)+^(^-^)] (-5) g-In^)+o(e), 

(78) 

Finally we consider the four-gluon vertex. Again we have a leading colour contribution and a 
contribution proportional to Nf. Representative diagrams are shown in fig. (11 e) and fig. (11 f). 
The ultraviolet subtraction terms are given by 
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Integration leads to 



Cr.nf = ^(2rt V °-^-^V p )^(^-ln^)+0(e). (80) 
3.4 The insertion operator L 

The integrated versions of the subtraction terms for the integrand of a one-loop amplitude have 
to be added back. The sum of all these subtraction terms is combined with the ultraviolet coun- 
terterm Sl^y and defines the insertion operator L: 

*<°>*L*<°> = 2Re^»)*(4'» + ^» + ^> + ^'»). (81) 
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In massless QCD this operator is rather simple and given by 
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Combining the insertion operator L with the insertion operator I we obtain in the massless case 
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We see that the sum of the insertion operators I and L is free of poles in the dimensional regular- 
isation parameter e. 

Let us now consider the massive case. It is sufficient to discuss the case of QCD amplitudes 
with one heavy flavour, the generalisation to several heavy flavours is straightforward. There are 
a few modifications. We have to take into account the heavy quark field renormalisation constant, 
which is given by 



Z 2 , e = i + ^c F ^~-4 + 31n^+o(oJ). 



(84) 



Secondly, the mass of the heavy quark is renormalised. For the heavy quark mass we have to 
choose a renormalisation scheme. In the on-shell scheme the mass renormalisation constant is 
given by 
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In the MS-scheme the mass renormalisation constant is simply given by 



(86) 



The renormalisation of the mass and of the heavy quark field will modify the ultraviolet coun- 
terterm and the insertion operator L. In order to present the counterterm related to the mass 
renormalistion it is convenient to define the quantity ®(°)(/?i, ...,p n ,g 7 m) through 
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We consider a QCD amplitude with n external partons, out of which n g are gluons, n q are mass- 
less quarks, nq are massless anti-quarks, hq are massive quarks with mass m and Hq are massive 
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anti-quarks with mass m. Obviously we have n q = nq and hq = hq We find for the operator L 
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Here we defined for heavy quarks 

Ye = Yg = o. 

The coefficient Cm depends on the renormalisation scheme for the mass and is given by 
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4 Contour deformation 

In this section we present the algorithm for the contour deformation. The guiding principles 
of the method are outlined in subsection 4.1. The individual parts are subsequently discussed in 
detail: Subsection 4.2 describes the details of the deformation. In order to resolve the direction by 
which we approach singular points we introduce additional Feynman parameters. Issues related 
to the numerical stability of the Monte Carlo integration are discussed in subsection 4.3. 



4.1 Overview of the contour deformation 

Having a complete list of ultraviolet and infrared subtraction terms at hand, we can ensure that 
the integration over the loop momentum gives a finite result and can therefore be performed in 
four dimensions. However, this does not yet imply that we can simply or safely integrate each of 
the four components of the loop momentum k M from minus infinity to plus infinity along the real 
axis. There is still the possibility that some of the loop propagators go on-shell for real values 
of the loop momentum. If the contour is not pinched this is harmless, as we may escape into 
the complex plane in a direction indicated by Feynman's +/8-prescription. However, it implies 
that the integration should be done over a region of real dimension four in the complex space 
C 4 . There are many choices for the contour deformation which are formally correct, in the sense 
that they avoid the poles of the propagators where a deformation is possible and deform the 
contour in the right direction. However, not all choices are suitable for a numerical Monte Carlo 
integration. Most choices will lead to large cancellations between different integration regions 
and therefore to large Monte Carlo integration errors. Finding an algorithm which determines a 
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correct deformation of the integration contour and which leads to a numerical stable integration 
is a challenging problem. We use several steps and techniques to achieve this goal. The critical 
regions in the integration domain are the regions where one or more propagators go on-shell. 

1. In a first step we introduce Feynman parameters as additional integration variables. This 
transforms the singular variety defined by the vanishing of some of the n propagators into 
a simple variety consisting of a single cone. For a single cone the required deformation for 
the loop momentum k can easily be stated. The deformation vanishes at the origin of the 
cone. 

2. In a second step we deform the Feynman parameters into the complex plane. This is 
necessary because the deformation of the loop momentum vanishes at the origin of the 
cone. The deformation of the Feynman parameters vanishes whenever Landau's equations 
are satisfied. 

3. The total deformation vanishes therefore if the loop momentum k is at the origin and if the 
Feynman parameters satisfy Landau's equations. Let us now consider an integrable singu- 
larity corresponding to the vanishing of n s i ng propagators in an amplitude with n external 
particles. The Feynman parametrisation has the disadvantage that it raises the denominator 
to the power n. This leads to numerical instabilities in the Monte Carlo integration for large 
values of n. In order to improve the stability of the Monte Carlo integration we expand the 
true propagators around propagators where we have added a small mass /j^ R . We will take 
/uj R purely imaginary with Im fj^ R < 0. (Comparing /uj R with jj^y we first note that we take 
them both purely imaginary. In practise we choose Im /jy V < Im /jf R < 0.) 

We remark that the singular region defined by the vanishing of the propagator of the UV- 
counterterm does not introduce additional problems: By choosing /jy V along the negative imag- 
inary axis we can ensure that this singular region is moved away from the integration contour. 

Let us now discuss how the contour deformation is done in detail. We would like to evaluate 
numerically the expression 

A d) _ A W _ A (i) _ A (i) _ A (!) - [J°L(gW _ G (1) -G (1) -G W ) (91) 

^subtr — ^bare ^soft ^coll ^UV _ J (2%) D V bare soft co11 uv / ' ^ ' 

The subtraction terms G^L and G^y are local counterterms. Therefore the final result is 
finite and the integral can be performed in 4 dimensions. We can write the integrand in the form 




P(k) and P\jv(k) are polynomials in the loop momentum k. Gi' and G I contribute only to the 
first term of the r.h.s of eq. (92), while g|j V contributes only to the second term on the r.h.s of 
eq. (92). On the other hand, G^L contributes to both terms on the r.h.s of eq. (92). The number 
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nuv is a positive integer, which by construction of the ultraviolet counterterms is not larger than 
four. The difference on the right-hand side of eq. (92) falls off at least like l/\k\ 5 for large \k\, 
which corresponds to the fact that the difference is UV-finite. The two individual terms 



P{k) 



Puv(k) 
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(93) 



This corresponds to the fact that all contributions are 



fall off at least like l/\k\ 2 for large \k\ 
maximally quadratic divergent. 

Eq. (92) shows clearly all possible singularities of the integrand. The integrand has singu- 
larities when one of the propagators \/(kj — m 2 ) goes on-shell as well as when the propagator 

of the ultraviolet counterterm l/(k — A*uv) § oes on-shell. The singularities of the ultraviolet 
counterterm are all localised on a single cone k — /j^y = and rather easy to handle. With the 
definition of P(k) and 7°uv(£) we can write the integral as 
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(94) 



In the next steps we define an integration contour which avoids the poles at k 2 — m 2 = and 
P-A/2 V = 0. 



4.2 Feynman parametrisation 

Let us now consider an integral of the form 
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(95) 



where P(k) and Z\jv(£) are polynomials in the loop momentum k, R(k) is a rational function in 
the loop momentum. The integration is over a complex contour in order to avoid - whenever 
possible - the poles of the propagators. The direction of the deformation is indicated by the 
+z'8-prescription. We now proceed as follows: We first introduce Schwinger parameters. We 
then deform the loop momentum and the Schwinger parameters into the complex plane. Each 
deformation will introduce a Jacobian. We can then integrate out one Schwinger parameter 
and we arrive at a formula equivalent to the Feynman parametrisation. We could have started 
directly with the Feynman parametrisation. The detour through the Schwinger parameters has 
the advantage that we obtain the correct Jacobian in a simple way. 
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Assume that A is a complex number with a positive imaginary part. Then we have 

oo 

-i j dte iAt . (96) 



1 

A 



We use this identity to rewrite the propagators with the help of the Schwinger parametrisation: 
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(-?')" j dt\... J dt n exp I i £ tj (kj - 

oo V j=1 



(97) 



For the argument of the exponential function we have 
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(98) 



We now set 



1 A 



(99) 
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We note that k as a function of the variables tj is homogeneous of degree 0. The Jacobian is 



-Ai. 



(100) 



Then we have 



n j n n 



(101) 



(=ij=i 



where denotes the Euclidean scalar product. The imaginary part of eq. (101) vanishes for 
k = 0. The resulting singularities can still be avoided by deforming the Schwinger parameters 
into the complex plane. Therefor we set 



tj + iltjfij (fx,..., i n ) , (3 j ft,..., ?„) 



k=\ 



n I n 
a=l \J>=1 



(102) 



such that all tj are real and positive. We note that tj as a function of the variables t\, t n is 
homogeneous of degree 1 . We denote by J the Jacobian of this transformation 



J(t\,...,t n 



dti 



(103) 
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Note that the Jacobian is a homogeneous function of degree in the variables t\, t n . Since 
R(k) and ...,?„) are both homogeneous of degree in the variables t\, t n we can integrate 
out one variable. We make use of the identity 



J d"tf(t) = f d n t f dth(t-t) f(t)= f d"xd( 

tj>0 tj>0 t>0 xj>0 ^ 
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with Xj = tj/i. We also write 
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(106) 



This is the standard formula for the Feynman parametrisation supplemented with the correct 
Jacobian corresponding to the deformation given in eq. (102) and eq. (104). The deformation for 
the parameters xj vanishes for xj = 0. The deformation has the additional property that it does 
not necessarily vanish in regions where two parameters are equal. 

We further note that / is a function homogeneous of degree (— n) in the variables x\, 
x n . We may use this property to replace the integration over the (n — 1) -dimensional simplex 
by an integration over the ^-dimensional hyper-cube. This is discussed in appendix B.2. The 
integration over the loop momentum W is over four real variables from minus infinity to plus 
infinity. It is convenient for a numerical Monte Carlo integration to map this region onto the 
finite region [0, l] 4 . This is discussed in appendix B.l. 

It remains to discuss the ultraviolet subtraction terms. For the arbitrary four- vector occur- 
ring in the ultraviolet subtraction terms we can make the choice 
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Since the Feynman parameters x; are deformed into the complex space it follows that Q 1 is in 
general a four- vector with complex entries. With this choice of Q M the relation between and 
W is given by 



1 



(108) 
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For the propagators of the ultraviolet subtraction terms we have then 

k 2 -nly = Ukok-nly. (109) 

The Euclidean norm ko k is always non-negative. If we set ^ v purely imaginary with Im /j^y < 
we ensure that this quantity is never zero. 

4.3 Improving the numerical stability 

In the integrand of eq. (105) we have the function 

j n n 

L(k;xi,...,x n ) = 2ixkok+ — Y^Y, XaSabXb - ( 110 ) 

lX a=\b=\ 

Let us discuss the conditions under which this function vanishes. To this aim we focus on the 
imaginary part. In order to simplify the following discussion we multiply this function by x. The 
factor x is never zero and therefore uncritical. We find that to order X the imaginary part of (xL) 
is given by 



n / n \ 

L xj I Sjixi 

7=1 ' V=l / 
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n / n 

L L Sabh 
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The expression on the right-hand side is always non-negative. The expression is zero if k = 
and if for all j 6 {1 , . . . , n} we have either 



xj = (112) 



or 

itsjixi = 0. (113) 

i=i 

Eq. (112) and eq. (1 13) are the well-known equations of Landau. These two equations together 
with the condition k = give therefore necessary conditions for the vanishing of the function 
L. Let us now consider a subset S C {1, ...,«} of rc s j ng elements, such that eq. (112) is fulfilled 
for all j G {l,...,n}\5 and eq. (113) is fulfilled for all j G S. For n s i ng < n this corresponds to 
a non-leading Landau singularity. In the final formula for the contour deformation in eq. (105) 
the function L occurs to the power (— n), but only a power of (— n s j ng ) can be attributed to the 
underlying physical configuration. The remaining (n — n smg ) powers are artificially introduced 
by the Feynman parametrisation. These powers are compensated by the integration over the 
Feynman parameters in the directions corresponding to {1, ...,n}\5. However this integration is 
done numerically with Monte Carlo methods. The net effects are numerical instabilities for large 
n from the regions where the function L is close to zero. 
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In order to improve the efficiency of the Monte Carlo integration we have either to avoid that 
the function L gets close to zero or to reduce the power n to which the function L is raised. We 
present here an efficient method which is based on the first strategy. In appendix C we describe 
a method based on the second strategy. 

In order to avoid that the function L gets close to zero we start from the identity 

T(n m + n) ( V' 1R 
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r ( n ) mf^o r ("IR + 1 ) V L - XflfR ) 



(114) 



Here we introduced an additional parameter /jj R , which we will take purely imaginary with 
Im /j^ R < 0. This ensures that 

Im (L-jc/Jj R ) > 0. 

It follows that the expression (L — x/^ R ) 1S never zero - 

The sum on the right-hand side of eq. (1 14) converges for 
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< 1. 



To leading order in X we can show that we always have 



< 1, 



(116) 



(117) 



and that the equal sign implies L = 0. For L = the sum on the right-hand side of eq. (114) is 
divergent, but in this case at the same time L 11 on the left-hand side is an ill-defined expression. 
We then truncate the series at order /Vir. This amounts to the replacement 
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(118) 



The net effect is that the critical expression L n is replaced by (L — ^j R ) _ " _ " IR , where (L—x/u^) 
is always non-zero. The replacement depends on two parameters /jj R and Nir, which we can use 
to control the quality of the approximation. We note that the replacement in eq. (1 18) is exact in 
the limit fjj R — > 0. 



5 Additional remarks 

In this section we offer a few additional remarks, which might help to avoid possible pitfalls. 
Subsection 5.1 is devoted to diagrams like massless tadpoles, which give zero in an analytical 
calculation. For this reason they are often discarded in an analytical calculation from the very 
beginning. Discarding these diagrams in a numerical calculation will spoil the local cancellation 
of the singularities, so these diagrams have to be kept. Subsection 5.2 discusses the scheme- 
independence of the infrared subtraction terms. Subsection 5.3 is devoted to subleading tree-level 
partial amplitudes, which can occur in the infrared subtraction terms. 
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Figure 12: Diagram (a) shows the top-level diagram for the left-moving primitive qggq ampli- 
tude. Diagrams (b) and (c) show lower-point diagrams contributing to the same amplitude. In an 
analytical calculation within dimensional regularisation these diagrams yield zero after integra- 
tion. In a numerical calculation these diagrams have to be included in order not spoil the locality 
of the subtraction terms. 



5.1 Locality of the subtraction terms and vanishing diagrams 

A one-loop amplitude can be represented by Feynman diagrams. The Feynman diagrams con- 
tributing to an amplitude are the ones contributing to the corresponding amputated Green's func- 
tion. For a one-loop amplitude all Feynman diagrams which are not self-energy insertions on 
external lines contribute. Fig. (12) shows some diagrams contributing to the left-moving prim- 
itive qggq amplitude. We focus here on the diagrams (b) and (c) of fig. (12). We first note that 
diagram (b) is not a self-energy insertion on the external line with momentum p2- This diagram 
belongs to the amputated Green's function. In an analytical calculation using dimensional reg- 
ularisation the diagrams (b) and (c) yield zero after integration. For this reason they are often 
discarded in an analytical calculation from the very beginning. This cannot be done in a numeri- 
cal calculation. Diagram (b) yields zero because of an exact cancellation between an ultraviolet 
divergence and an infrared divergence. Neglecting this diagram in a numerical calculation would 
spoil the local cancellations of the singularities. 

Diagram (c) yields zero because the massless tadpole corresponds to a scaleless integral. 
The corresponding integral is set to zero within dimensional regularisation. Diagram (c) has 
a quadratically divergent ultraviolet singularity. In a numerical calculation this singularity is 
subtracted out with the ultraviolet subtraction terms for the gluon propagator. 

5.2 Scheme independence 

In this subsection we investigate the scheme independence of the infrared subtraction terms. 
Contrary to the subtraction terms for the real emission, there is no dependence on the variant of 
dimensional regularisation (conventional dimensional regularisation, 't Hooft-Veltman scheme, 
four-dimensional scheme) in the integrated result for the infrared subtraction terms for the loop 
integrand. At first sight this is puzzling: In the real emission case the variant of dimensional 
regularisation introduces a scheme-dependent finite term. On the other hand unitarity requires 
that this scheme-dependent finite term cancels in the final result. It is instructive to investigate 
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how this cancellation occurs. For simplicity we discuss the massless case. The solution comes 
from the LSZ reduction formula: The renormalised one-loop amplitude with n g gluons, n q quarks 
and nq antiquaries is related to the bare amplitude by 

Aren(Ph-,Pn,a s ) = {z^^ (z]'^ A hwt (p h ...,p n ,Zp^ VV.) ■ (119) 
Z g is the renormalisation constant for the strong coupling, given by 

z * - 1+ s(-t)I +0 «)- (120) 

Z2 is the quark field renormalisation constant and Z3 is the gluon field renormalisation constant. 
The LSZ reduction formula instructs us to take as field renormalisation constants the residue of 
the propagators at the pole. In dimensional regularisation this residue is 1 for massless particles 
and therefore the field renormalisation constants are often omitted from eq. (1 19). However Z2 = 
Z3 = 1 is due to a cancellation between ultraviolet and infrared divergences [11]. In Feynman 
gauge we have 

Z 2 = l + ^C F (- L\ +0 (a$), 

A% \e IR e uv J 

Z 3 = l + ^-(2C A -%)(^---L)+o(a 2 s ). (121) 

471 \Eir Euv J 

Here we indicated explicitly the origin of the 1/e-poles. These poles introduce scheme-dependent 
finite terms of ultraviolet and infrared origin with opposite sign. Now the cancellation of the 
scheme-dependent finite term of infrared origin is as follows: The scheme-dependent finite term 
from the real emission contribution cancels with the scheme-dependent finite term of infrared 
origin from the renormalisation constants. This leaves the scheme-dependent finite terms of 
ultraviolet origin in the renormalisation constants and in the bare one-loop amplitude. These 
remain and give the result of the calculation in the chosen scheme. One can convert from one 
scheme to another by a finite renormalisation. We remark that the scheme dependence of the 
bare one-loop amplitude is entirely of ultraviolet origin. 

5.3 Subleading tree-level partial amplitudes 

In this subsection we discuss the tree-level partial amplitudes, which occur in the infrared sub- 
traction terms. We would like to point out, that a tree-level partial amplitude occurring in the 
infrared subtraction terms need not be identical to the ones occurring in the colour decomposition 
of the leading order calculation. This is best explained through an example. Let us consider the 
left-moving one-loop primitive amplitude A^ l \qi,g2,q3, £4) with the cyclic order gi,g2,<?3,g4- 
The top-level one-loop diagram for this amplitude is shown in the left diagram of fig. (13). 
This primitive one-loop amplitude has soft singularities, when either the gluon in the loop with 
momentum k\ or the one with momentum ^2 becomes soft. The amplitude has collinear singu- 
larities when either £4] \k\, k\\ \k% or k%\ I&3. In all these cases the infrared subtraction terms are 
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Figure 13: The left diagram shows the top-level diagram for the one-loop primitive qgqg ampli- 
tude. The tree-amplitude entering the infrared subtraction terms is given by the sum of the two 
diagrams on the right. 




Figure 14: The diagrams contributing to the partial amplitude A (°\q i,g2 ,£4,<?3) ■ 



proportional to the tree-level partial amplitude with the cyclic order <7i,g2,<73 1( ?4- The two tree 
diagrams contributing to this tree-level partial amplitude are shown in the middle diagram and in 
the right diagram of fig. (13). Note that these diagrams only involve the quark-gluon vertex, but 
no three-gluon vertex. On the other hand the colour-decomposition of the leading order tree-level 
amplitude is given in eq. (27) and involves the partial amplitudes with the cyclic order 

A (0) (qi,g 2 ,84,<h) , ^ (0) (<Zi,£4,£2,? 3 ) • (122) 

In these two partial amplitudes the quark and the anti-quark are always adjacent in the cyclic 
order. Each of these partial amplitudes consists of two Feynman diagram. For the amplitude 
^■ (#1 » £2j £4 » 93) me diagrams are shown in fig. (14). The corresponding diagrams for the 
amplitude A (°>(q 1,^4, g2, #3) are obtained by exchanging g2 ^ £4- 

The tree-level partial amplitude A*- ) (qi,g2, <?3, £4) can be calculated (as any other tree-level 
partial amplitude) with the recursive techniques of subsection 2.5. However, if the partial am- 
plitudes in eq. (122) are already known, the partial amplitude A^°\qi,g2,q3, gA) is simply given 
as 

A^(q hg2 ,q3,g4) = A(%i,g 2 ,£4,43)+A (0) (<7i,g 4 ,£2,<73). (123) 
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In the sum on the right-hand side the diagrams with the non-abelian three-gluon vertex drop out, 
due to the anti-symmetry of the colour ordered three-gluon vertex. 

This is easily generalised to amplitudes with more gluons: For an amplitude with the cyclic 
order q,gi,...,g m ,q,g m+h ...,g n we have 

A (0) (^,gl,...,gm,?,gm+l,---,gn) = £ A (0) (q,g ai , ...,ga„,q) ■ (124) 

ae{l,...,m}UI{n,...,m+l} 

This expresses the amplitude with the cyclic order q,g\, ...,g m ,q,g m +i,...,g n in terms of ampli- 
tudes with the cyclic order q,g ai , in which the quark and the anti-quark are adjacent. 
The sum is over all shuffles of the set {1, ...,m} with the set {n, ...,m+ 1}. A shuffle is a permu- 
tation of the n elements, which preserves the relative order of the elements in the first set, as well 
as the relative order of the elements in the second set. In order to prove the relation in eq. (124) 
one observes that in the sum on the right-hand side all non-abelian vertices involving gluons 
from both sets drop out. Let P a and Pf, denote momenta, which are the sums of momenta of a 
subset from {p\,...,p m }. Similar, we denote by Pj a momentum, which is the sum of momenta 
of a subset from {p m +\, •••>/?«}• In the sum over shuffles each three-point vertex connected to 
lines with the (ordered) momenta P a ,Pi will also occur in the order Pi,P a . Due to the antisym- 
metric nature of the colour ordered three-point vertex this will add up to zero. For the four-point 
vertices we have a similar situation. Four-point vertices with mixed momenta occur always in 
the combination P a ,Pb,Pi, in the combination P a ,Pi,Pb as well as in the combination Pj,P a ,Pb. 
Again, the three combinations will add up to zero. 

6 Checks and examples 

In this section we discuss a few checks and examples. One of the simplest processes where our 
method can be applied are the NLO corrections to the process 7* — > qq. This tests the locality of 
the subtraction terms and is discussed in subsection 6.1. In subsection 6.2 we describe a strong 
check for the contour deformation. We study the accuracy by which we can evaluate numerically 
the one-loop three point function with massless internal propagators and massive external lines. 
This integral is finite and we compute this integral for a contour based on a process with n 
external legs. We can think of this integral as occurring in an amplitude with n external legs 
where (n — 3) propagators are pinched. 

6.1 The process 7* — > qq 

The NLO corrections to the process 7* — > qq is a simple example where our method can be 
applied. We discuss this process in detail. This example illustrates nicely how the subtraction 
terms cancel the divergences locally. 

The Born amplitude squared for 7* — > q(p\) + q(p2) is given in D = 4 — 2e dimensions by 

= 4N c Q 2 q e 2 (l-E)s, (125) 
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with 5 = 512. Here we use the notation sy = (/?,- + Pj) 2 - The elementary electric charge is denoted 
by e and Q q gives the electric charge of the quarks in units of the elementary electric charge e. 
The virtual correction is given with k\ = k — p\, hi = k — p\ — p2 and £3 = k by 
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Using momentum conservation and the Dirac algebra in D dimensions this is equivalent to 
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In an analytical calculation all terms in the last line will yield zero after integration. The infrared 
subtraction terms are given by 
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As before we have set k = k— Q. As ultraviolet subtraction term we need only a subtraction term 
for the quark-photon vertex. The form of this subtraction term is identical to the subtraction term 
for the subleading colour contribution for the quark-gluon vertex. The latter has been given in 
eq. (74). We have 



,(0)%(i) 

l 3 -^UV 
2 



(129) 



(0) 



C F g 2 2R&S~ l p 2z 



d D k 



S12 



(2 Pl -k) {2 P 2-k)-2zk 2 -^l w 



The subtracted one-loop amplitude is finite and we may take the limit D — > 4. We obtain 
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This integral is finite and within the numerical approach it will be evaluated with Monte Carlo 
techniques. For the example discussed here, the integral is rather simple and can also easily be 



40 



evaluated analytically. The analytical result for this integral is given by 
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We can therefore compare the numerical result with the exact analytical result. With the methods 
of section 4 we easily achieve with a numerical Monte Carlo integration a precision at the per 
mille level for the integral in eq. (130). 

The subtraction terms for the integrand of the one-loop amplitude are added back in integrated 
form. For the process y* — > qq we have in addition that there are no ultraviolet counterterms from 
renormalisation, i.e. 
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We therefore have 
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This defines the insertion operator L for the process y* qq as 
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Let us also consider the real correction. The Born amplitude squared for y* — > q{p\) +g(p2) + 



(0) 



c F ^v c e^v8(i-£) 



^-2^-2^ + (l-£)^ + (l-£)^-2£ 
^12*23 S\2 *23 *23 *12 

\2 



where we used the notation sn^ = (pt + pj + pk) . The dipole subtraction terms read 
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For the real corrections one evaluates numerically the finite integral 



(0) 



£>qg,q ®qg,q 



SC F N c Q 2 q e 2 g 2 J d<h 



Sl2 



S23 



SU+S13 ^13 + 523 



■?12(S13+S23, 

(136) 
(137) 



+ 8 
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in four dimensions. The unresolved phase space measure is given in D dimensions by 

(4tz) £ - 2 



r(l-e) 



(*123) 



1-e 



J d 3 xd ^1 — Y* x i \ x \ 



— £ — £ — £ 

x 2. X 3 ' 



(138) 



where 



512 *23 *13 
XI = , X 2 = , X3 



S123 



*123 



*123 



(139) 



For the simple example considered here the integral in eq. (137) can also be done analytically. 
The result reads with s = 5123 



(0) 



®qg,q 'Dqg.q 



(0) 



\ 1 



(140) 



^^unres 

Integrating the dipole subtraction terms analytically in D dimensions one finds 

S £ l /J 2£ f J^unres {®qg,q + ®qg,q) = (141) 

+ 0(6). 



(0) 



1 



£ 2 £ 



+ - 3-2hw + ln z ^-3hw + 10--7^ 



We may now sum up the integrated subtraction terms from the virtual part and the real part. This 
yields the expression 

2Re 4 0) * (4¥ + ^soft + *S + *$) + Se V 2£ / ^unres ( -Dqg,q + <Dqg,q) = 



2% 



C F 



(0) 



10 -3 Re In ) + O (e) , (142) 
j"uv. 



which corresponds to the sum of the insertion operators I + L. We thus have for the process 
f->qq 



I + L 



2k' 



C F 10 - 3 Re In 



A'uv 



+ £ 



(143) 



Note that this expression is finite. 

In summary the complete NLO corrections to the process 7* — > qq is the sum of three contri- 
butions, which are given in eq. (130), eq. (137) and eq. (142), respectively. Each contribution is 
individually finite. Summing them up one obtains the well-known result 



2Re^ 0) % (1) +S-V 2e / 



(0) 



3 a, 
2 ' 27i 



C F 



.1 



(0) 



+ o (e) . (144) 
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Figure 15: The diagram on the left shows the top-level diagram for a process with n external 
legs. Pinching (n — 3) loop propagators results in a three-point function, shown on the right. 



6.2 Precision and accuracy of the contour integration 

The previous example shows nicely the local cancellation of all divergent parts. However, the 
virtual correction of the previous example has only three external legs and the resulting numer- 
ical integration in eq. (130) is rather trivial. Complications due to a large number of external 
particles do not enter the previous example. In order to test our method in situations with a large 
number of external particles we start by considering a configuration with n external particles. 
The top-level one-loop diagram associated with this configuration is a diagram with n internal 
loop propagators. If we now pinch (n — 3) of the n internal loop propagators, we obtain a one- 
loop three-point function. This is illustrated in fig. (15). The three-point function has the external 
momenta 

Pl=pi + ...+Pi-l, P 2 = pi + ...+pj- h P 3 = Pj + ...+p n . (145) 
We now consider the massless one-loop scalar three-point function 

9 r d 4 k 1 

1 = 1671 / T^a-. 5 9 (146) 

J (2%) 4 i k 2(k-P l ) 2 (k-P 1 -P 2 ) 2 

in the case P\ ^ 0, P 2 ^ and P 2 ^ 0. This integral is finite and the analytical result is well 
known [117-119]. With the notation 

hl= p2_p2_pl 82= p2_p2_p2 h= p2_ p 2_p2 

A 3 = (Pf) 2 + (Pj) 2 + (Pj) 2 - 2P\pI - 2P\P\ - 2P]P\ , ( 1 47) 



the three-mass triangle / is given in the region P\,P\,P 2 < and A3 < by 
2 



X 



(148) 



Cl 2 ( 2arctan ( ^y^H +C1 2 (^arctan ( ^-^1 +C1 2 |^2arctan ( 
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CI2 (x) denotes the Clausen function. In the region Pf.P^P 2 < and A3 > as well as in the 
region P\,P\ < 0, P\ > (for which A3 is always positive) the integral / is given by 



Re 



2(Li 2 (-pJt) +Li 2 (-py)) +ln(pjc)ln(py) +ln ( J - ) 



m0(/f) 



In 



1 + px 
1+PJ 



(61 + V^T) (83 + V^3) 



A 3 I (8l - y/A^) (83 - V 7 ^) 



7T 2 ' 

+ y 



(149) 



where 



P 2 
'2 



"^9> ^=7^ P 



/> 2 



n2 
r 3 



2P 2 



(150) 



The step function Q(x) is defined as 0(x) = 1 for jc > and 0(jc) =0 otherwise. Equipped with 
the analytical result we may now test the precision and accuracy of the numerical integration. We 
can think of the three-point function as a finite lower-point function occurring in the calculation 
of an 77-point amplitude. Therefore we have to check if the numerical integration is precise 
and accurate with respect to a contour which avoids all the singular surfaces defined by the n 
propagators 



(k - qi 



0, 



l,..,n. 



(151) 



and not just the three propagators 



k 2 = 0, 



(fc_ J p 1 ) 2 = 0, (k-P l -P 2 ) 2 = 0. 



(152) 



The former is more challenging than the latter. In general, the three-point integral / evaluates to 
a complex number. We test the real part and the imaginary part separately. 

We start with a configuration corresponding to just n = 3 external particles and then increase n 
step by step. All numerical integrations are done with the help of the Vegas algorithm [120, 121]. 
After a warm-up phase with five iterations of 10 5 evaluations each, we obtain the numerical re- 
sult from 20 iterations of 10 evaluations each. Up to n = 5 we find no problems in obtaining 
the correct answer from the numerical Monte Carlo integration with a precision better than one 
per cent by using just the methods of section 4.2 alone. However, starting from n = 6 we find 
numerical instabilities (and as consequence large statistical errors) by using just the methods of 
section 4.2. The situation is significantly improved with the help of the method of section 4.3. 
This method introduces two parameters /j 2 r and Nir. We recall that we choose /j 2 R as a purely 
imaginary number with Im /j 2 r < 0. It is convenient to parametrise fj 2 R by a dimensionless num- 
ber T|ir through 



a4 



(153) 



where Q is the centre-of-mass energy of the process under consideration. Fig. (16) shows the 
result of the numerical integration normalised to the analytical result for n = 6 external momenta 
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Figure 16: Results for the numerical integration of the three-point function normalised to the 
analytical result for a contour based on six external particles for various values of T]ir and Mr. 
The errorbars indicate the statistical error of the Monte Carlo integration. 
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for various values of T|ir and Nir. We have plotted the results for T|ir = 0.04, T|ir = 0.06 and 
TliR = 0.08 as a function of Nir. In all cases the correct value is approached as we increase Afo. 
As expected, we observe that the series converges faster for smaller values of T|ir. However, as 
can be inferred from the errorbars in fig. (16) smaller values of T|ir imply larger statistical errors. 
By a suitable choice of T|ir and A^ir one obtains again a precision better than one per cent. We 
have tested the numerical accuracy and precision of the three-point function with a contour up to 
n = 9 propagators in the loop. 

7 Conclusions 

In this paper we have presented a complete algorithm for the numerical calculation of one-loop 
QCD amplitudes. The algorithm consists of two main ingredients. First a set of subtraction terms, 
which approximate locally the soft, collinear and ultraviolet divergences of one-loop amplitudes. 
The integrals over these subtraction terms are simple and calculated analytically in D dimensions 
once and for all. The difference between the loop amplitude and the subtraction terms is finite 
and can be evaluated numerically in four dimensions. 

Secondly, the algorithm consists of a method to deform the integration contour into the com- 
plex space. Here we add for a process with n external particles to the four dimensions of mo- 
mentum space n additional dimensions corresponding to the Feynman parameters. This leads to 
a simple and general prescription for the contour deformation. In order to improve the stability of 
the numerical Monte Carlo integration we expand all propagators in the loop around propagators 
where a small fictitious width has been added. 

The algorithm is formulated at the amplitude level and does not rely on Feynman graphs. 
Therefore all required ingredients can be calculated efficiently using recurrence relations. The 
algorithm applies to massless partons as well as to massive partons. 

Although we have presented here the algorithm for QCD amplitudes, the extension to the 
electroweak sector and the Higgs sector, as well as to physics beyond the standard model seems 
to be straightforward. 
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A Colour ordered Feynman rules 

In this appendix we give a list of the colour ordered Feynman rules. They are obtained from the 
standard Feynman rules by extracting from each formula the coupling constant and the colour 
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part. The propagators for quark, gluon and ghost particles are given by 

. + m 



k 2 — m z 



ig 



k 2 

i 



(154) 



The colour ordered Feynman rules for the vertices are 



-if, 



(155) 



ft* 



f v (k$-k\)+x^ 



v 3 ""2 




= l 



iki. 



(156) 



B Generating the random points 

The numerical integration discussed in section 4 is over the four-dimensional loop momentum 
space and an (n — 1) -dimensional Feynman parameter space. After the contour integration we 
integrate over real variables. The integration over the four-dimensional loop momentum space 
extends for each dimension from minus infinity to plus infinity. The integration region over the 
Feynman parameters is an (n — 1) -dimensional simplex embedded into an /^-dimensional space. 
In order to be able to use standard methods for the Monte Carlo integration we map these domains 
to the unit hyper-cube. 
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B.l Generating the loop momentum 

In this appendix we show how to rewrite the integral 

J d 4 kf(k), 



(157) 



where the integration is over four real variables Icq, k\, &2 an d £3 from minus infinity to plus 
infinity, as an integral over the finite domain [0, l] 4 . We use four uniformly in [0, 1] distributed 
random numbers u\, 114 and define the four quantities kg, ^, and § by the equations 



kg 



n 

tan -mi : 



(^ — sin^cos^) — U2 = 0, 



= arccos (1 — 2^3) 
<|> = 2nu 4 . (158) 

Hi is an arbitrary scale. The second equation is solved numerically for ^. We then set 



sin sin (|) 

^o = ^£cos^, fc r = &£sin^, k = k r I sin0cos(|) 



(159) 



COS0 



The (inverse) Jacobian of this transformation is 



p(k) = 



du 



dk 



7T- 



+k?y+v 



^0 



(160) 



We then have 



J d 4 kf(k) 



j d 4 u 

[0,1] 4 



P (k) 



(161) 



B.2 Generating the Feynman parameters 

In this appendix we consider an integral of the form 




(162) 



over the standard (n — 1) -dimensional simplex. We assume that the function f(x\, ...,x n ) is ho- 
mogeneous of degree (— n) in the variables xj. We may then replace the integration over the 
simplex by an integration over the n faces of the 77-dimensional hyper-cube which are not coor- 
dinate subspaces. These n faces are defined by 

xi = 1, <*, <!, i^l. (163) 
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We can rewrite the integral as 



i 11 
/ = £ fd n x8(xi-l)f{x u ...,x n ) = f t [d n x8( Xl -l)f( Xl ,...,x n ) fdlnl"- 1 . 



Setting 



ui 



(164) 



we arrive at 



I = n j d n uf(x\,...,x n ) 



(165) 



where 



Ui 



* 7 

J max 



(166) 



and j max is the index such that Uj mRX > Ui for all i. 



C Reducing the power in the denominator 

In this appendix we show how to reduce the number of propagators, and in consequence also 
the power to which the denominator function L occurs in eq. (105). We split the integration into 
different channels. In each channel the number of propagators is reduced. In different channels 
different propagators are removed. Let us consider a process with n pr0 cess external partons. We 
denote the number of loop propagators in a specific channel by n c hannel- F° r each channel we 
have again the kinematical n cha nnei x n channe i-matrix and the (n channe i - 1) X (n channel - 1)- 
dimensional Gram matrix G$ . These kinematical quantities are obtained from the kinematics 
for the process with n pr0 cess propagators by pinching (n pr0 cess - ^channel) propagators. 

The reduction technique exploits the fact that in four dimensions we can have maximally four 
linear independently vectors and is well known from analytical calculations. Here we use it in 
a purely numerical context. The technique makes use of a special property of the kinematical 
matrix S defined in eq. (36), namely that for each channel we can always find numbers b\, £2, 
b„ , , such that 

"channel 



^channel 



£ Sijbj = 1 for all 1 < i < n c hannel- 
7=1 



(167) 
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As a standalone method the reduction technique is not as efficient as the method presented in 
section 4.3, mainly because the coefficients b\ are in general large and alternating in sign. How- 
ever the reduction technique can be combined with the method of section 4.3 in a way we will 
indicate below. This improves the efficiency of the method of section 4.3 even further. 

There is a minor modification which we have to make to the method presented in section 4.2: 
The reduction technique introduces additional ultraviolet subtraction terms. These additional 
ultraviolet subtraction terms cancel when summed over all channels. However the cancellation 
requires that the four- vector Q M is chosen to be the same in all channels. The choice in eq. (107) 
does not meet this requirement, since this choice depends on the Feynman parameters and there- 
fore on the channel. An alternative choice is 

-i "process 

& = ~ E process.,- (168) 

"process j = j 

We use the choice in eq. (168) in connection with the reduction technique. Again we can check 
that the ultraviolet propagators never go on- shell with a suitable choice of /j^y F° r me propaga- 
tors of the ultraviolet subtraction terms we have now 

P-A^v = 2ifk+^Vjofk+^-^VoV-i^jy + 2~k-V + V 2 , (169) 



where 



J "channel } "process 

V = ~ 52 ^channel, i — ~ H ^process J- (170) 

■* i = \ "process j = \ 

Since the Euclidean norm is always non-negative we can ensure for infinitesimal X that the imag- 
inary part of this expression is always positive by setting purely imaginary with 

? 1 

Im^uv < --maxVoV. (171) 

Doing so, we ensure that the propagators of the ultraviolet subtraction terms never go on-shell. 

In practise we use the reduction technique in combination with the techniques presented in 
section 4. We rewrite an integral of the form 

d A k R(k) 



(27!) 



4 ^process / 

n fe-mj 

7=1 V 



(172) 



d 4 k R(k) f d 4 k 



+ 



r d*k 

J J2k) 



( 27 r;)4 "process , ^ J ( 2 J 

n { k j- m j-viR 



R{k) R(k) 



^process / \ '^process / 

n {kj-mj) n (fi-m)-^ 



(173) 



The first term is calculated directly with the help of the Feynman parametrisation. To the second 
term we apply first the reduction identities and then the methods of section 4.2 and 4.3. 
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C.l Reduction identities 



In this appendix we show how to reduce in a channel with n propagators the number of propa- 
gators from n to (n— 1). In the reduction identities we distinguish the cases n > 6 and n < 6. In 
the former case there is a complete reduction to (n — 1), while in the latter case an integral with n 
propagators remains. It may seem that in the case n < 6 not much is gained, but it can be shown 
that the remaining integral with n propagators is better behaved in the infrared. We start with the 
case n>6, the case n < 6 will be discussed later. 

We will make use of the property that we can always find numbers b\, Z?2, b n such that 



E Stjbj 

7=1 



1 for all 1< i < n. 



(174) 



The kinematical n x n-matrix S was defined in eq. (36). The kinematical matrix S and the num- 
bers bj are independent of the loop momentum k. We have collected useful information for the 
computation of the numbers bj in the appendix C.2. For n > 6 the numbers b\, b n have the 
additional properties that 



I bj = 0, £ bjqj = 0. 



(175) 



This leads to the identity 



1 = ^{(kj-mf) 

i=i 



k 2 ~^w+ £ 2k-(q r -Q) 



j= l v 



(^ 2 -^v)" 



(176) 



The second term inside the curly bracket is an ultraviolet subtraction term and compensates the 
quadratic and linear growth of kf — m? with large \k\ . The complete expression inside the curly 
bracket goes to a constant for large \k\. The sum of all the extra ultraviolet subtraction terms 
vanishes. If we start from an integral of the form 



/ 



d 4 k R(k) 



n 



m' 



(177) 



where R(k) is a rational function in the loop momentum with poles only at 



2 

"A'UV 



0. 



(178) 



we obtain the following reduction: 



" r <Fk 



Ri(k) 



n 



m 



(179) 
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with 



Ri(k) = biR(k) { 1 



r=\,r^i 



(180) 



Eq. (179) splits an integral with n propagators into n integrals with [n— 1) propagators each. 
This reduces the number of propagators which need to be taken into account for the contour 
deformation. The rational functions Rj(k) appearing in eq. (179) have again only poles at the 
location of the ultraviolet subtraction propagator. We can iterate this procedure until we arrive at 
channels with five propagators. 

Let us now consider the case n < 5. Again we can find numbers b\, b n , such that 



E s a b i 



1 for all 1< / < n. 



But in contrast to the previous case we now have in general 



The reduction identity reads now 

d 4 k R(k) 



I = 



+ 1 



(fk 



Ri(k) 



n 



k)- m ) 



/ ' ' i I../-' 

where Ri(k) is defined as in eq. (180). The rational function R(k) is given by 



R(k) = R(k) 



C + Cuv 



n(k)- m ) 



C and Cuv are given by 



2\ x i 



C = -Bk 2 + 1 £2k-(b m ) + Y, 1 £ b i^ 2 j- 2< li-qj-mj) 

i=l i=l. ;=1 

C uv = B{P-^)+2k-(-{n-l)BQ + p(B-br)qr). 



(181) 



(182) 



(183) 



(184) 



(185) 



In the formula for the coefficient C the variables xj (j = l,...,n) can be arbitrary. We can take 
them equal to the Feynman parameters. The quantities B and x are given by 



B=Y,bi, x=Y,*i- 

i=\ i=l 



(186) 
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It is advantageous to slightly rearrange the formula for the coefficient C. We first write the loop 
momentum W as 

kf = K» + -Y Xi cfi, (187) 

where the Xj's are the Feynman parameters. In terms of the shifted loop momentum we have 
for the coefficient C 




(188) 



Eq. (183) does not eliminate the n-point integrals. Instead this reduction identity provides an 
additional factor C in the numerator for the n-point integrals. The factor C vanishes, whenever the 
denominator after Feynman parametrisation and contour deformation vanishes. This improves 
the numerical stability. 



C.2 The coefficients bi 

In this appendix we collect useful information related to the fact that the vector (1, 1, 1) is al- 
ways in the range of the kinematical matrix S [68,71,72, 122-124]. We recall that the kinematical 
matrix S was defined in eq. (36) asanx rc-matrix by 

Sij = (qi-qf) -mj-mj. (189) 

S depends only on the external momenta and the internal masses, but not on the loop momentum. 
It is a highly non-trivial statement, that we can always determine n numbers b\, b n such that 

n 

£ Sijbj = 1 for all 1 < i < n. (190) 
For generic external momenta the matrix S is invertible for n < 6. In this case bi is given by 

bi = t ( 191 ) 

7=1 

For n > 1 we always have det S = and therefore S is not invertible. However the vector ( 1 , . . . , 1 ) 
is still in the range of the linear map defined by S. In this case the solution for the numbers bi is 
not unique. A possible set of coefficients bj can be obtained from the singular value decomposi- 
tion of the (n — 1) x (n — 1) Gram matrix Gy — 2(qi — q a )(qj — q a ) where the indices take the 
values i,j ^ a. The starting point is to express the kinematical matrix Sy in terms of the Gram 
matrix: 

Sij = -G^+V^+V^, V^ = { qi - qa f-ml (192) 
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The condition in eq. (190) translates into the two equations 

£ df)bj = Biv^-v^ 

j^(vj a) -vi a) )bj = l-25vj fl) , (193) 

n 

where we set B = £ bj. An equation of the form Gx = y has a solution if and only if y = GHy, 

j=i 

where H is the pseudo-inverse of G. The pseudo-inverse H of a symmetric matrix G is uniquely 
defined by the properties 

HGH = H, GHG = G, GH = HG. (194) 
In the case where y = GHy holds the solution to the equation Gx = y is given by 

x = Hy+(l-HG)u, (195) 

where u is an arbitrary vector. If G does not have maximal rank one can show that y = GHy 

implies y = 0. Applied to our case this in turn implies 5 = 0. Noting that vj"^ — = S a j we 
therefore have the equations 

I Gifbj = 0, 

7=1 

To solve these equations one first determines (n — 1) numbers (b\, ...,b a -i,b a +i, ■■■,b n ), which 
define a vector in the kernel of G. The second equation is then used to fix the value of b a . The 
matrix G has the singular value decomposition 

4 

Gij = £t/*w,(V r ) (197) 

k=\ 

where U and V are orthogonal [n — 1) x (n — 1) matrices. The kernel of G is spanned by the 
vectors V^, V ; y„_i). We can therefore set 



bi = 


V i5 . 

W s > l<a ' 


bi+i = 


V i5 


b a = 


n 



(198) 
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II 



The normalisation factor W5 we get from £ S a jbj = 1: 

.7=1 



W5 = Efay-SaJV/s+E^-S^Vijs- 



(199) 




Let us now consider the four-vector r = £ /3,g ; . 



For n > 6 we can show that 



<7a-r = 



3V 



(200) 



for all a,b G {l,...,n}. In particular this implies 



q a -r 



q n -r = 0. 



(201) 



Therefore it follow that r = 0. 

Let us summarise the situation: We distinguish three cases. 

1. n < 5. For generic external momenta the kinematical matrix S is invertible and the coef- 
ficients bi are unique. The sum of all coefficients B is non-zero and the four-vector r is 
non-zero. 

2. n = 6. For generic external momenta the kinematical matrix S is invertible and the coeffi- 
cients bi are unique. The sum of all coefficients B is zero and the four-vector r is zero. 

3. n>l. The kinematical matrix S is not invertible and the coefficients bi are not unique. The 
sum of all coefficients B is zero and the four-vector r is zero. 

At the end of this appendix we consider the special case where a constant quantity fij R is added 
to all masses squared: 



The kinematical matrix Sij is real, whereas the matrix may be complex. Suppose we are 
interested in the coefficients b' t , which satisfy 




The quantity /j^ R may be complex, however we assume that the masses m,- are real. As before we 
have the kinematical matrix Sij defined by 





1. 



(205) 
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These can be obtained without complex arithmetic. We first determine the real coefficients b\ 
defined by 

Hsijbj = 1 (206) 

.7=1 

and set then 



One easily verifies that the coefficients in eq. (207) satisfy eq. (205). In particular we have b\ = bi 
for n>6. 
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